require(tidyr)
require(dplyr)
require(stringr)
require(rgeos)
require(rgdal)

out_path = "D:/Research/global_refor_mac/model_runs/"

model_list = readRDS(paste0(out_path, "glob_model_run_stand.rds"))

fold_path = "D:/Research/global_refor_mac/results/rds_reg_files/"

# upload simulation regressions and dataframe for all three continents
cont_df <-  readRDS(paste0(fold_path, "all_cont_reg.rds"))
cont_df = cont_df[which(is.na(cont_df$carbon_baccini) == FALSE), ]

# calculate bau emissions 
bau_position_list = 1
dollar_50_position = 7
dollar_20_position = 4

# calculate the deforestation emissions
def_emission_bau = model_list[[bau_position_list]][[3]]
def_emission_bau = apply(def_emission_bau[,2:4], 1, sum) / (cont_df$shp_size * 100)

def_emissions_50dol = model_list[[dollar_50_position]][[3]]
def_emissions_50dol = apply(def_emissions_50dol[,2:4], 1, sum) / (cont_df$shp_size * 100)

ref_emission_bau = model_list[[bau_position_list]][[6]]
ref_emission_bau = apply(ref_emission_bau[,2:4], 1, sum) / (cont_df$shp_size * 100)

ref_emissions_50dol = model_list[[dollar_50_position]][[6]]
ref_emissions_50dol = apply(ref_emissions_50dol[,2:4], 1, sum) / (cont_df$shp_size * 100)


def_emissions_20dol = model_list[[dollar_20_position]][[3]]
def_emissions_20dol = apply(def_emissions_20dol[,2:4], 1, sum) / (cont_df$shp_size * 100)

ref_emissions_20dol = model_list[[dollar_20_position]][[6]]
ref_emissions_20dol = apply(ref_emissions_20dol[,2:4], 1, sum) / (cont_df$shp_size * 100)


### create the output dataset 
outdata = model_list[[bau_position_list]][[3]]

# create graph data frame
graph_df = data.frame("uni_ID" = outdata$uni_ID,
                      "country" = outdata$country,
                      "def_bau" = def_emission_bau,
                      "re_df_50d" = def_emission_bau - def_emissions_50dol,
                      "ref_bau" = ref_emission_bau,
                      "eh_rf_50d" = ref_emissions_50dol - ref_emission_bau,
                      "re_df_20d" = def_emission_bau - def_emissions_20dol,
                      "eh_rf_20d" = ref_emissions_50dol - ref_emission_bau,
                      stringsAsFactors = FALSE)

# read in all of the unique gridded areas
str_matrix = str_split(outdata$uni_ID, "_", simplify = TRUE)
uni_split_grids = unique(paste0(str_matrix[,2], "_", str_matrix[,3]))

for(i in 1:length(uni_split_grids)){
  
  # read in the shapefile 
  shp_in = readOGR("D:/Research/global_refor_mac/data/unit_analysis",
                   paste0(uni_split_grids[i], "_fish_cl"))
  shp_in@data$uni_ID = paste0(as.numeric(row.names(shp_in@data)), "_",
                              uni_split_grids[i])
  
  # merge the files 
  shp_data = left_join(shp_in@data, graph_df, by = "uni_ID")
  shp_in@data <- shp_data
  
  # write out the shapefile 
  writeOGR(shp_in, 
           "D:/Research/global_refor_mac/results/shp_graph", 
           paste0(uni_split_grids[i], "_out_shp"), 
           driver="ESRI Shapefile", 
           overwrite_layer=T)
  
  
}


### In this last step we append all shapefile and rasterize the shapefile
### We create a map of global emissions reduction under a given carbon price

# create list of all shapefiles in directory
list_shp <- list.files(paste0("D:/Research/global_refor_mac/results/shp_graph/"), pattern = "\\.shp$")

# create character string of first shapefile in list
shp1 <- paste0(paste0("D:/Research/global_refor_mac/results/shp_graph/"), list_shp[1])

# use GDAL to create merged shapefile with the first shapefile
system(paste0("ogr2ogr ", paste0("D:/Research/global_refor_mac/results/shp_out/merged_shp.shp "), shp1))

# the loop append all other shapefiles in directory to the merged shapefile
# this creates a shapefile that includes all the areas from Asia, Latin America, and Africa
for(i in 2:length(list_shp)){
  
  # creates character string of shapefile that will appended to merged shapefile
  shp2 <- paste0("D:/Research/global_refor_mac/results/shp_graph/", list_shp[i])
  
  # append shapefile to merged shapefile
  system(paste0("ogr2ogr -update -append ", 
                paste0("D:/Research/global_refor_mac/results/shp_out/merged_shp.shp "), 
                shp2, " -nln merged_shp"))
  
}

# rasterize the merged shapefile and create raster image of shapefile 
#system("gdal_rasterize -a cum_emh -tr 0.05 0.05 -a_nodata -999 F:/WFWN_working_paper/data/graphics/emission_chart/2015_2020/raster/merged_shp.shp F:/WFWN_working_paper/data/graphics/emission_chart/2015_2020/raster/emissions_reduct_crb.tif")
system(paste0("gdal_rasterize -a def_bau -tr 0.05 0.05 -a_nodata -999 ", 
              paste0("D:/Research/global_refor_mac/results/shp_out/merged_shp.shp "),
              paste0("D:/Research/global_refor_mac/results/shp_out/def_bau.tif ")))

# rasterize the merged shapefile and create raster image of shapefile 
#system("gdal_rasterize -a cum_emh -tr 0.05 0.05 -a_nodata -999 F:/WFWN_working_paper/data/graphics/emission_chart/2015_2020/raster/merged_shp.shp F:/WFWN_working_paper/data/graphics/emission_chart/2015_2020/raster/emissions_reduct_crb.tif")
system(paste0("gdal_rasterize -a eh_rf_50d -tr 0.05 0.05 -a_nodata -999 ", 
              paste0("D:/Research/global_refor_mac/results/shp_out/merged_shp.shp "),
              paste0("D:/Research/global_refor_mac/results/shp_out/enhanced_rem_ref_50dol.tif ")))

# rasterize the merged shapefile and create raster image of shapefile 
#system("gdal_rasterize -a cum_emh -tr 0.05 0.05 -a_nodata -999 F:/WFWN_working_paper/data/graphics/emission_chart/2015_2020/raster/merged_shp.shp F:/WFWN_working_paper/data/graphics/emission_chart/2015_2020/raster/emissions_reduct_crb.tif")
system(paste0("gdal_rasterize -a re_df_50d -tr 0.05 0.05 -a_nodata -999 ", 
              paste0("D:/Research/global_refor_mac/results/shp_out/merged_shp.shp "),
              paste0("D:/Research/global_refor_mac/results/shp_out/avoided_def_50dol.tif ")))

# rasterize the merged shapefile and create raster image of shapefile 
#system("gdal_rasterize -a cum_emh -tr 0.05 0.05 -a_nodata -999 F:/WFWN_working_paper/data/graphics/emission_chart/2015_2020/raster/merged_shp.shp F:/WFWN_working_paper/data/graphics/emission_chart/2015_2020/raster/emissions_reduct_crb.tif")
system(paste0("gdal_rasterize -a ref_bau -tr 0.05 0.05 -a_nodata -999 ", 
              paste0("D:/Research/global_refor_mac/results/shp_out/merged_shp.shp "),
              paste0("D:/Research/global_refor_mac/results/shp_out/ref_bau.tif ")))

# rasterize the merged shapefile and create raster image of shapefile 
#system("gdal_rasterize -a cum_emh -tr 0.05 0.05 -a_nodata -999 F:/WFWN_working_paper/data/graphics/emission_chart/2015_2020/raster/merged_shp.shp F:/WFWN_working_paper/data/graphics/emission_chart/2015_2020/raster/emissions_reduct_crb.tif")
system(paste0("gdal_rasterize -a eh_rf_20d -tr 0.05 0.05 -a_nodata -999 ", 
              paste0("D:/Research/global_refor_mac/results/shp_out/merged_shp.shp "),
              paste0("D:/Research/global_refor_mac/results/shp_out/enhanced_rem_ref_20dol.tif ")))

# rasterize the merged shapefile and create raster image of shapefile 
#system("gdal_rasterize -a cum_emh -tr 0.05 0.05 -a_nodata -999 F:/WFWN_working_paper/data/graphics/emission_chart/2015_2020/raster/merged_shp.shp F:/WFWN_working_paper/data/graphics/emission_chart/2015_2020/raster/emissions_reduct_crb.tif")
system(paste0("gdal_rasterize -a re_df_20d -tr 0.05 0.05 -a_nodata -999 ", 
              paste0("D:/Research/global_refor_mac/results/shp_out/merged_shp.shp "),
              paste0("D:/Research/global_refor_mac/results/shp_out/avoided_def_20dol.tif ")))


# crop the shapefiles
system('gdalwarp -cutline D:/Research/global_refor_mac/data/continent/three_continent_dis.shp -crop_to_cutline -dstnodata -9999 D:/Research/global_refor_mac/results/shp_out/ref_bau.tif D:/Research/global_refor_mac/results/shp_out/ref_bau_final2.tif') 
system('gdalwarp -cutline D:/Research/global_refor_mac/data/continent/three_continent_dis.shp -crop_to_cutline -dstalpha D:/Research/global_refor_mac/results/shp_out/avoided_def_50dol.tif D:/Research/global_refor_mac/results/shp_out/avoided_def_50dol_final.tif') 
system('gdalwarp -cutline D:/Research/global_refor_mac/data/continent/three_continent_dis.shp -crop_to_cutline -dstalpha D:/Research/global_refor_mac/results/shp_out/enhanced_rem_ref_50dol.tif D:/Research/global_refor_mac/results/shp_out/enhanced_rem_ref_50dol_final.tif') 
system('gdalwarp -cutline D:/Research/global_refor_mac/data/continent/three_continent_dis.shp -crop_to_cutline -dstalpha D:/Research/global_refor_mac/results/shp_out/def_bau.tif D:/Research/global_refor_mac/results/shp_out/def_bau_final.tif') 
system('gdalwarp -cutline D:/Research/global_refor_mac/data/continent/three_continent_dis.shp -crop_to_cutline -dstalpha D:/Research/global_refor_mac/results/shp_out/avoided_def_20dol.tif D:/Research/global_refor_mac/results/shp_out/avoided_def_20dol_final.tif') 
system('gdalwarp -cutline D:/Research/global_refor_mac/data/continent/three_continent_dis.shp -crop_to_cutline -dstalpha D:/Research/global_refor_mac/results/shp_out/enhanced_rem_ref_20dol.tif D:/Research/global_refor_mac/results/shp_out/enhanced_rem_ref_20dol_final.tif') 


# crop the shapefiles
system('gdalwarp -cutline D:/Research/global_refor_mac/data/continent/three_continent_biome_arabia.shp -crop_to_cutline -dstnodata -9999 D:/Research/global_refor_mac/results/shp_out/ref_bau.tif D:/Research/global_refor_mac/results/shp_out/ref_bau_final_b.tif') 
system('gdalwarp -cutline D:/Research/global_refor_mac/data/continent/three_continent_biome_arabia.shp -crop_to_cutline -dstnodata -9999 D:/Research/global_refor_mac/results/shp_out/avoided_def_50dol.tif D:/Research/global_refor_mac/results/shp_out/avoided_def_50dol_final_b.tif') 
system('gdalwarp -cutline D:/Research/global_refor_mac/data/continent/three_continent_biome_arabia.shp -crop_to_cutline -dstnodata -9999 D:/Research/global_refor_mac/results/shp_out/enhanced_rem_ref_50dol.tif D:/Research/global_refor_mac/results/shp_out/enhanced_rem_ref_50dol_final_b.tif') 
system('gdalwarp -cutline D:/Research/global_refor_mac/data/continent/three_continent_biome_arabia.shp -crop_to_cutline -dstnodata -9999 D:/Research/global_refor_mac/results/shp_out/def_bau.tif D:/Research/global_refor_mac/results/shp_out/def_bau_final_b.tif') 
system('gdalwarp -cutline D:/Research/global_refor_mac/data/continent/three_continent_biome_arabia.shp -crop_to_cutline -dstnodata -9999 D:/Research/global_refor_mac/results/shp_out/avoided_def_20dol.tif D:/Research/global_refor_mac/results/shp_out/avoided_def_20dol_final_b.tif') 
system('gdalwarp -cutline D:/Research/global_refor_mac/data/continent/three_continent_biome_arabia.shp -crop_to_cutline -dstnodata -9999 D:/Research/global_refor_mac/results/shp_out/enhanced_rem_ref_20dol.tif D:/Research/global_refor_mac/results/shp_out/enhanced_rem_ref_20dol_final_b.tif')
